function output = gen_data_NL(input_mat, pibar)
T = input_mat.T;
lag = input_mat.lag;
rho = input_mat.rho;
pi = input_mat.pi;
pi_res = input_mat.pi_res;
theta = input_mat.theta;
X = input_mat.X;
Y = input_mat.Y;
Z1 = input_mat.Z1;
sigma_v = input_mat.sigma_v;
sigma_u = input_mat.sigma_u;

% Generate the data
Ystar = zeros(T, 1);
Xstar = zeros(T, 1);
vstar = zeros(T, 1);
ustar = zeros(T, 1);

for t = 1:T
    if t <= 2
        Xstar(t) = X(t);
        Ystar(t) = Y(t);
    else
        xi1 = randn(1, 1);
        xi2 = randn(1, 1);
        vstar(t) = sigma_v * xi1;
        ustar(t) = sigma_u * (rho * xi1 + sqrt(1-rho^2) * xi2);
        Xstar(t) = [1, Z1(t, :), Ystar(t-2)] * pi + [1, Z1(t, :), Ystar(t-2), Z1(t, :).^2, Ystar(t-2)^2] * pi_res * pibar + vstar(t);
        Ystar(t) = [1, Xstar(t)] * theta + ustar(t);
    end
end

% Define Z2
if lag >= 2
    Z2 = zeros(T-lag, lag-1);
    for l = 2:lag
        Z2(:, l-1) = Ystar(1+lag-l:(end-l));
    end
    z = [Z1((lag+1):end,:) Z2];
else
    z = [Z1((lag+1):end,:)];
end
% 
y = Ystar((lag+1):end);
x = Xstar((lag+1):end);
% 
output.y = y;
output.x = x;
output.z = z;
end